Mechanism for bipolar resistive switching in transition metal oxides 
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We introduce a model that accounts for the bipolar resistive switching phenomenom observed in 
transition metal oxides. It qualitatively describes the electric field-enhanced migration of oxygen 
vacancies at the nano-scale. The numerical study of the model predicts that strong electric fields 
develop in the highly resistive dielectric-electrode interfaces, leading to a spatially inhomogeneous 
oxygen vacancies distribution and a concomitant resistive switching effect. The theoretical results 
qualitatively reproduce non-trivial resistance hysteresis experiments that we also report, providing 
key validation to our model. 



PACS numbers: 73.40.-c, 73.50. -h 
I. INTRODUCTION 

There is a great deal of experimental activity currently 
devoted to explore new technologies for the next gen- 
eration of electronic memory devices^. Among various 
promising options, the resistive random access memory 
(RRAM), which is based on the resistive switching (RS) 
phenomenon, has emerged as a preeminent candidate^. 
The RS effect is a large, reversible and nonvolatile change 
in the resistance after the application of voltage or cur- 
rent pulses. The typical RRAM system has a capacitor- 
like structure composed of insulating or semiconducting 
materials sandwiched between two metal electrodes. RS 
has been observed in a wide variety of systems, suc h as 
simple and complex oxides, organic compounds, etcP^. 
However, there are specific characteristics of the switch- 
ing behavior observed in each type of material. In the 
case of binary oxides, which are highly insulating it is 
believed that the RS effect may be due to the formation 
and rupture of conductive filaments within the insulating 
medicP^. In contrast, in the more conducting or semi- 
conducting complex oxides with perovskite structures, 
such as doped cuprates and manganites, the relevance 
of oxygen vacancies is often invoked. Despite a burst- 
ing body of experimental data that is rapidly becoming 
available^!!] th e precise mechanism behind the physical 
effect of RS remains elusive. A few qualitative models 
have been proposed emphasizing differen t aspects: elec- 
tric field-induced defect migratio n 1 5 * 13 * 14 !, phase separa- 
tion^, tunneling acros s int erfacial domains^, control of 
Shottky barrier's heighPEll, etc. A general consensus has 
emerged on the empirical relevance of three key features: 
(i) a highly spatially inhomogeneous conduction in the 
low resistive state, (ii) the existence of a significant num- 
ber of oxygen vacancy defects and (iii) a preeminent role 
played by the interfaces, namely, the regions of the oxide 
that are near each of the metallic electrodes which of- 
ten form Schottky barriers. Here, we shall introduce and 



study the behavior of a simple model that incorporates, 
at a qualitative level, those three features. By means 
of a numerical simulation we shall show that the model 
correctly reproduces key non trivial hysteresis cycles ob- 
served in experiments on perovskite- type transition metal 
oxides (TMO). 



II. MODEL 

Several experiments revealed that the conduction in 
the low resistive state is highly inhomogeneous and dom- 
inated by one dimensional paths that are associated with 
enhanced conduction channels^! These paths would 
be created upon an initial application of strong electric 
fields, that bring the dielectric close to its breakdown 
point. Thus, we shall assume that the electric trans- 
port is dominated by a single conductive path embedded 
within a more insulating host. 

The second important feature incorporated into our 
model is the relevance of defects within the dielectric. 
Several experim ents point to a preeminent role played by 
oxygen vacancie d 5 * 14 * 18 * 19 -*. Moreover, it is a universal and 
salient feature of TMO that their resistivity is dramati- 
cally affected by the precise oxygen stoichiometry. One 
may thus expect that oxygen vacancy concentration may 
be the most significant parameter controlling the local re- 
sistivity, p, of a given material. This feature is included 
in our model by assuming that each nano-domain of the 
path is characterized by a certain concentration of oxy- 
gen vacancy defects, S. We adopt the most simple linear 
relation p oc <5, which follows from the fact that in TMO 
perovskites the presence of oxygen vacancies severely dis- 
rupts the electronic conduction properties. Nevertheless, 
we emphasize that the specific form of p{5) is not crucial 
for the results that we shall describe later on. 

The third important feature that our model incorpo- 
rates is the key role played by the interfaces^, namely, 
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FIG. 1: a) Schematic model with a single conductive channel 
within the dielectric. The three regions L, R and B corre- 
spond to the two high resistance interfaces and the more con- 
ductive central bulk, respectively. The small boxes indicate 
the domains, b) Detailed scheme of the conductive path. The 
grayscale qualitatively depicts the variation in the oxygen va- 
cancy concentration through the channel (darker corresponds 
to higher concentration). The top figure shows the initial 
state with uniformly distributed oxygen vacancies, and the 
bottom one shows the inhomogeneous distribution after the 
first few "forming" voltage cycles (see text). 



the regions of the dielectric that are in physical proximity 
to the metallic electrodes. There is growing evi dence th at 
these are the regions where the RS takes plac d 13 * 1 ^^"^. 

Our model is schematically shown in Fig.[T]and consists 
of a single conductive channel within a more insulating 
dielectric, which is represented by a one dimensional re- 
sistive network on N links. The first and last N] links 
correspond to high resistance interfacial regions next to 
the external electrode, and the central N — 2Nj links de- 
scribe the bulk section. Each link is characterized by a 
certain concentration of oxygen vacancies, which deter- 
mines the resistivity of the link. They may be physically 
associated to small domains of nanoscopic dimensions 
which may actually correspond to grains of the polycrys- 
talline oxide. We take pi = A a Si, with a = B if i is in 
the bulk (Ni < i < N - iV», a = L if i is in the left 
interface (i < JVj) and a — R if i is in the right inter- 
face (TV - Nj < i < N). In our study we set N = 100 
and N] = 10. The following equation specifies how the 
vacancies diffuse through the network domains under an 
external voltage, 



Pab = S a {l - 5 b ) exp(-V + AV a ) 



(1) 



is proportional to the concentration of vacancies present 
in domain a and to the concentration of "available va- 
cancy sites" at the target domain. The Arrhenius factor 
exp(— Vq), is controlled by a dimensionless constant, Vq, 
related to the activation energy for vacancies diffusion. 
The important factor exp(AV a ) models the enhancement 
(or suppression) of the diffusive process due to the local 
electric field at domain a. 

From Eq.Q, a constant Vq leads to an initially con- 
stant distribution = 8 for all i. The value of the 
initial constant concentration, 8 must be much smaller 
than one, since it physically represents the concentration 
of defects (oxygen vacancies) within a domain. We adopt 

S = io- 4 . 

Similarly to actual resistive switching experiments, 
we simulate the applied voltage protocol V(t) by lin- 
ear ramps that follow the sequence — > +V max — > — > 
—Vjnax — > (our convention is that the right electrode 
is grounded). The duration is of s time steps. The se- 
quence may be repeated a number of cycles n, for a to- 
tal duration r max = ns. In our simulations we choose 
Vo/Vmax = 0.016, that provides a non-negligible but slow 
diffusive contribution to the evolution of Si with respect 
to the total time duration of the simulations (ie, the to- 
tal number of time steps). We set V max = 1000 that 
provides for a sufficiently large electric stress. Our quali- 
tative results are rather robust with respect to the choice 
of model parameters, a detailed systematic study of their 
dependence is left for future work. 

The coefficients As, An and Al still remain to be spec- 
ified. With no loss of generality, we fix the value of the 
bulk coefficient to unity Ab = 1 and leave the interfacial 
An and Al free. In this initial study we shall concen- 
trate in the symmetric case, Ar = Al, that corresponds 
to most common experimental devices. 

The numerical simulations are performed through the 
following steps: (i) at each simulation time step t (1 < 
t < T max ) a given external voltage V(t) is applied to the 
resistive network (the electrodes are assumed perfect con- 
ductors). The current through the system is computed 
as I(t) = V(t)/R.T, with Rt the total (ie, two terminal) 
resistance 
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(2) 



It gives the probability for transfer of vacancies from do- 
main a to a nearest neighbor domain b. The probability 



where a = R, B, L denotes the three network regions, 
and c denotes an unessential geometrical constant related 
to the dimensions of the domains (which we set to unity) . 
(ii) we compute the local voltage profile Vi(t) = I{t)pi 
and the voltages drops AVi(t) = Vi + \{t) — Vi{t). (iii) we 
use Eq. ([!]) to compute all the oxygen vacancy transfers 
between nearest neighboring domains, and update the 
values Si(t) to a new set of concentrations §i{t + 1). (iv) 
we use these new values to recompute the current at r t +i 
under the applied voltage Vit + 1) as indicated in the 
first step. 
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FIG. 2: Top panel: Resistive hysteresis loop (Rt vs V) for 
the subsequent voltage cycles number 5 , 6 and 7. (a). Inset: 
the Rt vs I hysteresis (loop number 5) shows qualitatively 
similar results (the vertical axis scale is the same as in the 
main figure). Bottom panel: experimental hysteresis loops 
measured in manganite (b) and cuprate (c) devices. The first 
one was pulsed in current control mode while the second in 
voltage mode. The details of the experimental procedures are 
described in Refs. 1211 and 25, respectively. 



III. NUMERICAL RESULTS AND 
COMPARISON TO EXPERIMENTS 



We now turn to the discussion of our results. We set 
= Al = 1000 >> Ab = 1, which is consistent with 
our experimental data and with previo us repor ts in bulk 
and thin films of conducting perovskite d 14 * 21 * 2 ^. In Fig. [2] 
we show the results for the hysteresis loop of the total 
resistance Rt- The different data curves in Fig. [2^l dis- 
play the results in subsequent voltage cycles 5, 6, and 
7. The inset of Fig. [2ji, show that the hysteresis as a 
function of current remains qualitatively similar. We 
note that during the first few initial cycles the resistance 
shows non-repetitive memory effects that converge to a 
hysteresis loop with a stable shape. Interestingly, this 
is reminiscent of the initial "forming" that experimen- 
tal samples seem to require in order to start displaying 
reproducible switching effects. The peculiar type of hys- 
teresis loop that we obtain has been already reported 
by the Houston grourP^in experiments on (Pr,Ca)Mn03 
manganite systems, where it has been termed "table with 
legs" . It is evidently a non-trivial effect and we have ex- 
perimentally reproduced it in both, a related manganite 
(Pr,La,Ca)MnC>3 and a cuprate (YBa2Cu307_2;) sample, 
as shown in the bottom panels of Fig. [2j 

There are several features worth pointing out: During 
each voltage protocol loop, there is a clear variation of 
the resistance Rt between a rather broad maximum for 




FIG. 3: Density concentration profiles normalized to the ini- 
tial uniform density value 5i/5 , at the beginning of voltage 
cycles number 5, 6 and 7 in the symmetric system. 



a large range of V (ie, the "table"), and two relatively 
narrow minima (ie, the "legs"). These maximum and 
minima correspond to the high and low resistance states, 
R^ 1 and Rk , The Rt{V) loops are approximately sym- 
metric in V which reflects the left-right symmetry of the 
system. Throughout the voltage loop, the system begins 
in the initial R^ 1 state and undergoes the sequence of re- 
sistance changes — > R 



R T !I+ under positive bias, 



and then — > R^ u ~ — > R^ 1 under negative bias. The 
final state, at zero bias is R^ 1 , very close but not iden- 
tical to the previous initial state. Interestingly, a similar 
small drift is also observed in the experimental data. 

The qualitative agreement between our prediction of 
the "table with legs" and the experimentally observed 
hysteresis loops provides a significant validation for our 
model. Therefore, to gain physical insight into the mech- 
anism of the RS effect we shall discuss in detail the evolu- 
tion of the vacancies distribution under the applied volt- 
age protocol. In Fig. [3] we show successive snapshots of 
the oxygen vacancy concentration profile Si at the be- 
ginning of loops 5, 6 and 7. Recalling that the initial 
equilibrium distribution of oxygen vacancy concentration 
is uniform Si = S a , these curves reveal that, under the 
action of the repetitive voltage cycling, the Si evolve to- 
wards a new stable distribution. The salient features of 
the profile Si are a significant depletion in the interfacial 
regions and a strong accumulation peaks at both internal 
boundaries between the bulk and the interfacial regions. 
The reason can be understood as follows. The largest 
electric fields occur at the two interfacial regions since 
initially their resistance is much larger than the bulk one 
(At, Ah >> Ab)- Therefore, oxygen vacancy migration 
is enhanced in the interfacial regions, with the ions mov- 
ing either towards the electrodes or the bulk, depending 
on the direction of the applied voltage. When the ions 
reach the metallic electrodes they start to pile-up^. On 
the other hand, the vacancies that migrate towards the 
bulk eventually leave the interfacial region and enter the 
bulk. There, their diffusion virtually stops, since the elec- 
tric fields in the more conducting bulk are much smaller. 
Successive initial cycles yield a cumulative effect, with a 
depletion of the interfacial regions (and some pile-up at 
the edges) and a concomitant accumulation at the bulk 
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FIG. 4: Top (semi-log scale): Snapshots of the density concen- 
tration profiles for the symmetric configuration, normalized to 
the initial uniform density value 5i/S , during the hysteresis 
cycle (first half). Bottom: The corresponding profiles of local 
resistance pi for the same snapshots. Inset: Position of the 
snapshots in the (first half) hysteresis loop. 



side of the interfacial/bulk boundary. Significantly this 
accumulation, as seen in Fig.[3j is quite large and narrow 
leading to a substantial increase in the local resistance. 
This feature will be a key to understand the origin of the 
legs of the hysteresis loop, which we shall consider next. 

In Fig. [4] we show snapshots, during the first half of 
the voltage cycle, of the vacancy concentration and lo- 
cal resistance profiles in the interfacial regions and their 
boundaries with the bulk. These are the active regions 
of the system where the electric-field-enhanced migration 
takes place. Notice that since the current I(t) is uniform 
along the conductive path, the local electric fields are di- 
rectly proportional to the local resistance. Let's start at 
V=Q, from the state at the beginning of a cycle. Interest- 
ingly, the pi profile indicates significant electric fields in 
both interfacial regions which extends into the neighbor- 
ing bulk. Thus the depletion of vacancies in the inter- 
faces and the accumulation at the boundaries are such 
that compensate for the difference between the respec- 
tive A coefficients, to yield similar electric fields across 
the boundaries of the regions. Yet, at the very left end of 
the system there is a small pile-up of vacancies which, as 
soon as the voltage is ramped up, will translate into the 
largest local fields and initiate the ionic migration. At 
V/V mSiX =0.25, we observe a snapshot of the migration 
of the vacancies across the left interfacial region towards 
the bulk. From the resistance expression Eq. ([2| , so long 
the vacancies remain within the interfacial region, the 
system remains in 1 state. At V/V max =0.62, the va- 
cancies have moved out the left interfacial region and 
entered the bulk, where their migration suddenly stops. 



Once in the bulk, the contribution to the total resistance 
of these migrating vacancies is reduced (Ab « Al), 
thus the system reaches the R^P state (leg of the table) . 
Notice that in this state, the largest fields occur at the 
boundaries of the bulk and the interfacial regions. Thus, 
as V is further increased, the left interfacial region de- 
pletes further and therefore the voltage drop gets lower 
there. In contrast, on the right side the electric fields are 
further enhanced and there is now a migration of vacan- 
cies from the accumulation peak of the bulk towards the 
right interfacial region. This leads to an increase of the 
total resistance and, at the maximal voltage V/V max =l, 
we find that the vacancies have entered the interfacial 
region and already piled-up at the right end. Thus, the 
system is back to R^ 1 . The voltage protocol continues 
with the decrease of V back to zero but keeping the same 
(positive) polarity. Therefore no significant change in the 
Si and pi profiles occurs, and half of the table with legs is 
already formed. When the negative polarity part of the 
cycle begins a similar analysis follows, since the distribu- 
tion of vacancy concentrations is a mirror image of the 
initial one. This forms the other half of the table. 

To complete our study we consider the effect of intro- 
ducing an asymmetry in the model parameters Ar and 
Al. In Fig. [5] we display the results for the hystere- 
sis loops upon increasing the asymmetry. The results of 
the simulations show the gradual evolution of the "table 
with legs" towards a conventional rectangular hysteresis 
cycled Asymmetry in experimental samples can either 
be due to the fabrication process, as for instance, in the 
obvious case of using different type of metal for the elec- 
trodes. Though our samples were fabricated in a sym- 
metric configuration, we induced a significant asymme- 
try by vigorous pulsing with a given polarity (ie, poling) 
before measuring the resistance hysteresis characteristic. 
Significantly, as shown in Fig. [5j we find that the ex- 
perimental data obtained in an asymmetric manganite 
sample and the cuprate, induced by intensive same po- 
larity pulsing, are in good qualitative agreement with our 
simulations. These results provides additional validation 
to our model. 



IV. CONCLUSIONS 



To conclude, our results put on solid theoretical 
grounds the key role played by oxygen vacancies in the 
mechanism of resistive switching in TMO. They also pro- 
vide valuable insights, predicting a non-trivial spatial 
profile of the oxygen vacancy distribution which may be 
of help for device design. An exciting idea for future 
work is to explore the possibility of using atomistic, first 
principles, calculations to study properties of electrode - 
transition metal oxide interfaces to estimate the param- 
eters of the model and provide guidance in the material 
choice for actual memory devices. 
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FIG. 5: Resistive hysteresis loops obtained for an increasing degree of asymmetry. Al ~ 1000 and Ar = 1000, 100, 50 and 25 
in panels a, b, c and d respectively. The right panels show experimental data for a manganite (e) and a cuprate sample (f) 
that were rendered asymetric by means of intense and fixed polarity pulsing. The former was pulsed in current control mode 
and the latter in voltage control, similarly as in Refs.|2TJand 25 , respectively. Notice the qualitative similarity of data of panels 
e and f with panel d; and of data in Fig. 2b with panel a. The data of Fig. 2c seem to be in between the results of panel a 
and b. 
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